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Abstract. We considered a higher-dimensional extension for the replica-exchange Wang- 
Landau algorithm to perform a random walk in the energy and magnetization space of the 
two-dimensional Ising model. This hybrid scheme combines the advantages of Wang-Landau 
and Replica-Exchange algorithms, and the one-dimensional version of this approach has been 
shown to be very efficient and to scale well, up to several thousands of computing cores. This 
approach allows us to split the parameter space of the system to be simulated into several pieces 
and still perform a random walk over the entire parameter range, ensuring the ergodicity of the 
simulation. Previous work, in which a similar scheme of parallel simulation was implemented 
without using replica exchange and with a different way to combine the result from the pieces, 
led to discontinuities in the final density of states over the entire range of parameters. From our 
simulations, it appears that the replica-exchange Wang-Landau algorithm is able to overcome 
this difficulty, allowing exploration of higher parameter phase space by keeping track of the joint 
density of states. 


1. Introduction 

The Wang-Landau (WL) Monte Carlo method [Ij is a sampling technique that can be used to 
obtain the density of states (DOS) of a system. It has the property of generating a flat histogram 
in some random walk space, where the parameters for the random walk and the flatness criterion 
can be chosen according to the system of interest mm- Due to these characteristics and to the 
simplicity of the algorithm, the Wang-Landau approach has been studied and applied successfully 
to different kinds of systems (see, for examples, HI El El 13), including those with rough free 
energy landscapes and strong first-order phase transitions, which cannot be easily studied with 
traditional methods. 

Recently an efficient parallel scheme of the Wang-Landau method was implemented mm 
□BE], whereby independent random walks are carried out in different overlapping windows 
of energy and a replica-exchange scheme mm is used to exchange system configurations of 
neighboring energy windows. This parallel replica exchange Wang-Landau (REWL) method was 


shown to scale extremely well up to several thousands of computing cores. Therefore, it can be 
used to study more complex systems and much larger system sizes than it is possible with the 
serial Wang-Landau method. 

While the density of states as a function of a single parameter, such as the total energy, 
suffices to provide insight into many systems, there are a number of physical properties and 
systems whose investigation relies on a joint density of states dependent on additional parameters 
na cna em- Therefore, it would be of great interest to extend the one-dimensional REWL 
method [BJ EJ COI] to a higher-dimensional parameter space, and implement a similar method that 
combines enhanced sampling techniques liana na and the advantages of parallel computing 
to obtain a joint density of states. 

In this exploratory study we consider the Ising model, for which the density of states as a 
function of total energy is known exactly [15]. The phase space of this model, as a function 
of energy and magnetization, has cusps and empty regions with no realizable configurations. 
Simulations in this higher-dimensional parameter space are challenging because, depending on 
the sampling method used, the lower and higher energy edges of the phase space are hard to 
sample and thus have poor statistics. Moreover, the empty region between cusps can be hard 
to overcome in a simulation. 


2. Two-dimensional Replica-Exchange Wang-Landau sampling 

The REWL method [8] is arguably the most successful parallelization attempt to the WL method 
so far, and it consists of only three simple modifications to the serial version of the algorithm. 
Taking the Ising model as an example, a serial WL sampling of the total energy ( E ) and 
magnetization (M) phase space employs a two-dimensional random walk in the entire E and M 
ranges. In the process a random walker keeps track of the joint density of states g(E, M ) and 
the histogram of visited states H(E,M). Let A'(E',M') and A"{E",M") denote, respectively, 
the current and trial spin configurations of the system. The WL probability that the trial move 
is accepted is given by 


P(A' -> A") 


min 


' g(E',M') - 

. ’ g(E",M")_ 


(!) 


If A"(E", M") is accepted, the density of states g(E ", M") is adjusted by the modification factor 
/ > 1, that is, g(E",M") g(E",M") * /; and the histogram is updated with H(E",M") —>■ 
H(E ", M") + 1. If the trial state is not accepted, we update g(E',M') —> g(E',M') * / and 
H (. EM') —>■ H (. EM') +1. One Monte Carlo step is defined as Lx L trial state considerations. 

The random walker proceeds with choosing a new trial spin configuration and checking 
whether it is accepted, until H(E,M ) is considered flat, in accordance with a chosen criterion. 
At this point the current iteration is ended and the histogram is reset to zero. The next 
iteration begins with a finer /, reduced by some predefined rule (e.g. / —>• v7 ). More 

iterations are performed until the modification factor reaches a chosen final value /final; such 
as /final = exp(10 -8 ). The density of states is never reset and, at the end of the simulation, 
g(E, M) is expected to converge to the true value with an accuracy proportional to \//final- 

A replica-exchange Wang-Landau sampling starts by splitting the phase space into smaller 
overlapping subspaces (windows), as schematically shown in Fig. |T| For simplicity, only the 
energy range is split and the subspaces were artificially displaced vertically (in FigiTJb)) for 
clarity. 

Inside each subspace one or more independent random walkers perform a standard WL 
sampling, with independent g(E,M ) and H(E,M). After a predetermined number of Monte 
Carlo steps, a configurational exchange between two walkers in adjacent windows is proposed. 
The walker pairs are chosen at random. If either of them is not in the overlapping region, the 
exchange of configurations will not be attempted. However, if the walker pairs are in the same 







Figure 1. (color online) (a) Complete phase space for the Ising model on a 8 x 8 square lattice, 
(b) An example of splitting of the parameter space into overlapped, equal-sized energy windows. 
Energy windows are shifted upward progressively for clarity. 


overlapping region, the configuration swap between them (walkers i and j ) will be accepted with 
the probability: 


g r (E(X),M(X)) gj (E(Y),M(Y)) 

gi(E(Y),M(Y)) 9 j (E(X),M(X)) 


( 2 ) 


where X and Y denote, respectively, the spin configurations of walkers i and j before the 
exchange. Note that P acc satisfies detailed balance. All walkers (in all regions) move to the 
next iteration at the same time and the parallel simulation ends when the modification factor 
for every walker reaches /final- 

At this stage of the simulation, each walker is left with a joint DOS fragment (Figj2](a)) and 
one needs to combine all the fragments to recover g(E, M) over the entire phase space. To do so, 
we first sum over the magnetization and obtained the density of states as a function of energy: 


g(E)=J2 9 (E,M ), 

M 


( 3 ) 


as shown in Figj2](b). After that we calculate the derivative of ln[g(E)] as a function of E in 
the overlapping region. In this procedure, one finds a point E ma t where the difference between 
the inverse of the microcanonical temperature /3(E) = d \n[g{E)\/AE is a minimum. This E ma t 
point is used to match and adjust the two adjacent surfaces. That is, the whole DOS surface 
g(E, M) has to be rescaled at g(E mat , M) and the joint DOS is obtained by: 


g(E,M) 


gi(E,M), if E < .Emat, 
gj(E,M), if E > E mat . 


( 4 ) 


where i and j are, respectively, walkers in the lower and higher energy range and, gj(E,M ) is 
the shifted value of g(E, M ) for the jth walker. 

In principle, other ways of matching the 2D-DOS could be applied unum and other ways to 
split the phase space might be desired to increase the performance of the simulation ( see, e.g., 
Ref. [8j ). For example, one might be interested in splitting the magnetization instead of the 
energy range, or splitting both energy and magnetization ranges. In this case both dimensions 
of the overlapping phase space will need to be analyzed, one at a time, to find the E mat and 


































Figure 2. (color online) (a) Final g(E, M ) pieces of the two-dimensional Ising model on a 
8x8 square lattice. DOS pieces are shifted upward progressively for visualization purpose, (b) 
Individual pieces of g(E ) obtained by a summation over the magnetization range (Eq. ([3])). 
Vertical lines indicate the matching points (F mat ) for joining the DOS pieces together. F mat is 
determined as the point with minimal difference in the inverse microcanonical temperatures of 
the two DOS pieces. See text for more details. 


M mat matching points. Moreover, since multiple walkers can be used to sample a window 
simultaneously, it is possible to employ jackknifing or bootstrapping techniques to compute the 
statistical errors from a unique production run. 


3. Preliminary results for the two-dimensional Ising model 

From our simulations for the two-dimensional Ising model we examine (i) the generalization of 
the REWL approach to higher dimensions and ( ii ) the accuracy of the algorithm. We carried 
out serial 2D-WL and 2D-REWL simulations by performing random walks in total energy E and 
magnetization M space until the modification factor reaches /final = exp(10~ 8 ). We adopted the 
80% flatness criterion, i. e., all entries in the histogram H(E,M) are no less than the chosen 
percentage of the average height of the histogram. 

For the 2D-REWL simulations we divided the parameter space into equal-sized windows 
with restricted values of E and unrestricted M. We used a constant 50% overlap between 
neighboring windows and multiple random walkers to keep track of the joint DOS within each 
window. Figj3](a) shows the resulting DOS that is obtained in a few hours by 9 parallel walkers 
(the energy is split into 3 windows with three walkers per window). An equivalent serial WL 
sampling over the entire two-dimensional parameter space with a single random walker takes 
several days to complete. 

To compare the DOS obtained from WL and REWL sampling (denoted by g-wh(E,M) and 
5rf.wl (E, M), respectively), we calculated the relative difference: 


A di s(E,M) 


|ffw L (E,M) ~ 9rewl{E,M)\ 
g(E, M) 


( 5 ) 


where g(E , M ) = (gwi/E, M) + < 7 rewl(E', M))/2. The relative difference is plotted in Fig. EKb). 
Over the entire phase space, there is no apparent sign of systematic errors due to the matching 
of windows. The order of magnitude of the relative difference is comparable to the relative errors 
(as shown in Fig. IU(a)); they both share the same behavior of having larger values in the lower 
and upper energy edges. This is conceivable because of the rareness of configurations as well as 
the presence of cusps and empty regions in the phase space. 
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Figure 3. (color online) Logarithm of the complete joint DOS for the Ising model on a 16 x 16 
square lattice obtained by the two-dimensional REWL approach, (b) Relative difference between 
the joint DOS obtained from WL and REWL sampling. 


To verify the accuracy of the algorithm, we compute g(E ) with Eq. (0 so that we can compare 
with the exact results. FigU^a) shows the logarithmic DOS obtained by using 3 walkers per 
window in a 2D-REWL simulation along with the exact values, the relative errors are shown in 
the insert. We can see that the results obtained from our simulations agree extremely well with 
the analytic solutions. Figj4](b) shows the relative errors for the simulations employing 2 walkers 
per window and 20 walkers per window respectively. We observe that increasing the number of 
walkers, which update the DOS independently, decreases the sampling errors correspondingly. 




Figure 4. (color online) (a) Comparison of the density of states from 2D-REWL sampling and 
the exact results for a 16 x 16 square lattice (relative errors are shown in the insert), (b) Relative 
errors for a 8 x 8 lattice square using 2 walkers per window (upper panel) and 20 walkers per 
window (lower panel) respectively. 


4. Summary 

We implemented a two-dimensional replica-exchange Wang-Landau sampling method and 
applied it to determine the joint density of states g(E, M) for the Ising model on a square 
lattice. The joint DOS over the entire two-dimensional parameter space, obtained by combining 
the joint DOS in the different sampled windows, agrees very well with the results obtained with 






















a serial WL sampling over the entire two-dimensional parameter space. No artifacts in the final 
joint DOS were observed due to the parallelization and DOS joining procedures. We also showed 
that errors can be reduced by employing multiple random walkers in each window. Similar to 
the original REWL method, the higher-dimensional implementation of the REWL method is 
very efficient and it scales well with the number of CPUs. Hence, it can be used to study much 
larger systems than what is feasible using a serial WL method. 
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